clc;
clear all;
close all;

%% Load input image for edge detection algorithms
imgRGB = imread('retina4.jpg');
img = rgb2gray(imgRGB);
[rows cols] = size(img);

% default_edge_img = edge(img,'log');
% figure;imshow(default_edge_img);
% 
% default_edge_img = im2bw(default_edge_img);
% default_edge_img = remove_isolated_point(default_edge_img);
% default_edge_img = bwmorph(default_edge_img,'thin');

%% Laplacian of Gaussian
%% Adapt sigma
sigmas = [0.1 0.5 1 1.5 2 2.5 3 5 10];
hsizes = ceil(3 * sigmas);
thresh = 0.4;
figure;
% handle = tight_subplot(3,3,[.07 .01],[.1 .01],[.1 .01]);
for i = 1:length(sigmas)
    log_edge_img = ecen5283_edge_detect_log(img, sigmas(i), hsizes(i), thresh);
    subplot(3,3,i);
%     axes(handle(i));
    imshow(log_edge_img);
    title([' sigma =', num2str(sigmas(i)), ' ; size = ', num2str(2*hsizes(i)+1),' ; threshold = ', num2str(thresh)]);
end

%% Adapt threshold
sigmas = [2 2.5];
hsizes = ceil(3 * sigmas);
thresh = [0.3 0.35 0.4];
figure;
for i = 1:length(sigmas)
    for j = 1:length(thresh)
        log_edge_img = ecen5283_edge_detect_log(img, sigmas(i), hsizes(i), thresh(j));
        subplot(2,3,(i-1) * 3 + j);
        imshow(log_edge_img);
        title([' sigma =', num2str(sigmas(i)), ' ; size = ', num2str(2*hsizes(i)+1),' ; threshold = ', num2str(thresh(j))]);
    end
    
    
end

%% The best
sigmas = [1.9];
hsizes = ceil(3 * sigmas);
thresh = [0.3];
figure;
handle = tight_subplot(3,2,[.05 .01],[.03 .01],[.01 .01]); axes(handle(1));imshow(imgRGB);title('Original color image');
axes(handle(2));imshow(img);title('Gray image');
for i = 1:length(sigmas)
    for j = 1:length(thresh)
        [log_edge_img log_filtered_img zrc_img] = ecen5283_edge_detect_log(img, sigmas(i), hsizes(i), thresh(j));
        axes(handle(3)); imshow(log_filtered_img,[]);title('Filtered image by LoG filter');
        axes(handle(4));imshow(zrc_img,[]);title('Zero-crossing image');
        axes(handle(5));imshow(log_edge_img);
        title([' sigma =', num2str(sigmas(i)), ' ; size = ', num2str(2*hsizes(i)+1),' ; threshold = ', num2str(thresh(j))]);    
        %superimposed
        axes(handle(6));
        log_rgb = repmat(uint8(255*log_edge_img),[1 1 3]);        
        si_img = imadd(imgRGB, log_rgb);
        imshow(si_img);
    end        
end